快速傅里叶变换

一.复数

1.概念

复数就是形如 $a+bi$ 的数，其中 $a,b$ 是实数，且$b≠0,i^2=- 1$。

其中实数 $a$ 和 $bi$ 分别被称为复数的实部和虚部。

2.四则运算

1.加法 $(a+bi)+(c+di)=(a+c)+(b+d)i$

2.减法 $(a+bi)-(c+di)=(a-c)+(b-d)i$

3.乘法 $(a+bi) \times (c+di)=(ac-bd)+(ad+bc)i$

4.除法 $(a+bi) \div (c+di)=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{ac-bd}{c^2-d^2}+\frac{bc-ad}{c^2-d^2}i$

struct Complex {
    double x , y;

    Complex operator + ( const Complex &a ) const {
        return { x + a.x , y + a.y };
    }
    Complex operator - ( const Complex &a ) const {
        return { x - a.x , y - a.y };
    }
    Complex operator * ( const Complex &a ) const {
        return { x * a.x - y * a.y , x * a.y + y * a.x };
    }
    Complex operator / ( const Complex &a ) const {
        return { ( x * a.x - y * a.y ) / ( a.x * a.x - a.y * a.y ) , ( y * a.x - x * a.y ) / ( a.x * a.x - a.y * a.y ) };
    }
};

3.几何意义

每一个复数可以由直角坐标系上的点唯一确定，$x$轴称为实轴，$y$轴称为虚轴。

复数乘法有个重要的性质：模长相乘，幅角相加。

1.模长相乘

$l_1=\sqrt{a^2+b^2}$

$l_2=\sqrt{c^2+d^2}$

$l_3=\sqrt{(ac-bd)^2+(ad+bc)^2}$

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ =\sqrt{(ac)^2+(bd)^2-2abcd+(ad)^2+(bc)^2+2abcd}$

$~~~~~~~~~~~~~~ =\sqrt{(ac)^2+(bd)^2+(ad)^2+(bc)^2}$

$~~~~~~~~ =\sqrt{a^2(c^2+d^2)+b^2(c^2+d^2)}$

$=\sqrt{(a^2+b^2)(c^2+d^2)}$

$~~~~~~~~~~~=\sqrt{(a^2+b^2)}\sqrt{(c^2+d^2)}=l_1l_2$

2.幅角相加

证明 $\triangle AOP$ 与 $\triangle COB$ 相似即可。

首先有：

$\frac{AO}{CO}=\frac{OP}{OB}=\frac{1}{l_2}$

只需证：

$\frac{AP}{CB}=\frac{1}{l_2}$

$\frac{AP}{CB}=\frac{\sqrt{(a-1)^2+b^2}}{\sqrt{(ac-bd-c)^2+(ad+bc-d)^2}}$

$=\frac{\sqrt{a^2+b^2-2a+1}}{\sqrt{((ac)^2+(bd)^2+c^2-2abcd-2ac^2+2bcd)+((ad)^2+(bc)^2+d^2+2abcd-2ad^2-2bcd)}}$

$=\frac{\sqrt{a^2+b^2-2a+1}}{\sqrt{(ac)^2+(bd)^2+(ad)^2+(bc)^2+c^2+d^2-2ac^2-2ad^2}}$

$=\frac{\sqrt{a^2+b^2-2a+1}}{\sqrt{(c^2+d^2)(a^2+b^2-2a+1)}}$

$=\frac{1}{\sqrt{c^2+d^2}}=\frac{1}{l_2}$

二.复数单位根

1.概念

$n$ 次单位根是 $n$ 次幂为 $1$ 的复数。

$1$的模长为$1$，所以 $n$ 次单位根的模长也为$1$，即在单位圆上。

因为模长均为$1$，所以在单位圆上，我们只需考虑幅角。

不难想到，幅角为 $\frac{2\pi}{n}\alpha$ 的复数，都是单位根。

$n$ 次单位根一共有 $n$ 个，且等分单位圆。如图便是所有的 $6$ 次单位根。

类似的，$n$ 次单位根记为 $\omega_n^0$ ，$\omega_n^1$，$\omega_n^2$....$\omega_n^{n-1}$（逆时针标号）

2.性质

1.$\omega_n^0=1$

2.$\omega_n^k=(\omega_n^1)^k$

3.$\omega_n^i \times \omega_n^j=\omega_n^{i+j}$

4.$\omega_n^k=\omega_n^{k\%n}$

5.$\omega_{pn}^{pk}=\omega_n^k$

6.$\omega_n^{k+\frac{n}{2}}=-\omega_n^k$ ($n$ 为偶数)

7.$\omega_n^k=(\cos \frac{2k\pi}{n},\sin \frac{2k\pi}{n})$

8.$\sum_{i=0}^{n-1}\omega_n^i=0$

$\sum_{i=0}^{n-1}\omega_n^i=\frac{\omega_n^0(1-\omega_n^{n-1})}{1-\omega_n^1}=0$

三.多项式的表示

系数表示法

一个 $n-1$ 次多项式一共有 $n$ 项，每项前面的数字称为系数。

$f(x)=f_0+f_1x+f_2x^2+...+f_{n-1}x^{n-1}$

是最常见的一种表示方法，但乘法需要 $\Theta(n^2)$。

点值表示法

还记得周考卷的已知两点求直线解析式吗？

其实直线可以看做一个 $1$ 次多项式，我们知道两点即可得到解析式。

那么对于 $n$ 次多项式，只需要知道 $n+1$ 个取值即可唯一确定。

两个点值表示多项式的乘法只需要 $\Theta(n)$。

fft理论

既然点值表示法的乘法复杂度低，那么我们就可以将两个多项式由系数表示转化为点值表示，相乘后再还原为系数表示。

将多项式由系数表示到点值表示叫做离散傅里叶变换（$\text{dft}$），由点值表示到系数表示叫做逆离散傅里叶变换（$\text{idft}$）。

四.（逆）离散傅里叶变换（DFT/IDFT）

以下的 $n$ 均为 $2$ 的幂。

1.离散傅里叶变换

考虑一个 $n$ 项多项式 $f(x)$。

$f(x)=f_0+f_1x+f_2x^2+...+f_{n-1}x^{n-1}$

按奇偶分成两部分：

$f(x)=(f_0+f_2x^2....+f_{n-2}x^{n-2})+(f_1+f_3x^3+...+f_{n-1}x^{n-1})$

又有两个 $n/2$ 项多项式 $fl,fr$。

$fl(x)=f_0+f_2x+...+f_{n-2}x^{n/2-1}$

$fl(x)=f_1+f_3x+...+f_{n-1}x^{n/2-1}$

那么有: $f(x)=fl(x^2)+xfr(x^2)$

对于 $k(k < \frac{n}{2})$ 分别代入 $\omega_n^k$ ， $\omega_n^{k+n/2}$ 得到。

$f(\omega_n^k)=fl((\omega_n^k)^2)+\omega_n^kfr((\omega_n^k)^2)$

$f(\omega_n^k)=fl(\omega_{n/2}^k)+\omega_n^kfr(\omega_{n/2}^k)$

$f(\omega_n^{k+n/2})=fl((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}fr((\omega_n^{k+n/2})^2)$

$f(\omega_n^{k+n/2})=fl((\omega_n^{k+n/2})^2)-\omega_n^kfr((\omega_n^{k+n/2})^2)$

$f(\omega_n^{k+n/2})=fl(\omega_n^{2k+n})-\omega_n^k fr(\omega_n^{2k+n})$

$f(\omega_n^{k+n/2})=fl(\omega_n^{2k})-\omega_n^k fr(\omega_n^{2k})$

$f(\omega_n^{k+n/2})=fl(\omega_{n/2}^k)-\omega_n^k fr(\omega_{n/2}^k)$

然后我们会惊奇的发现，两个式子只有一个符号的区别。

现在就可以进行分治了。

2.逆离散傅里叶变换

不妨记 $G[k]$ 为多项式 $F$ 在 $\omega_n^k$ 的点值表示 ,则有

$G[k]=\sum_{i=0}^{n-1}F[i](\omega_n^k)^i$

那么代入 $G[i]$

$\sum_{i=0}^{n-1} (\omega_n^{-k})^i G[i]$

$\sum_{i=0}^{n-1} (\omega_n^{-k})^i \sum_{j=0}^{n-1}(\omega_n^i)^j F[j]$

$\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} (\omega_n^{-k})^i (\omega_n^i)^j F[j]$

$\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} \omega_n^{i(j-k)} F[j]$

令 $p=j-k$。

1.当 $p=0$ 时，即 $j=k$

$sum=\sum_{i=0}^{n-1}\omega_n^0F[k]=nF[k]$

2.当 $p\not=0$ 时

$sum=\sum_{i=0}^{n-1}\omega_n^{ip}F[k+p]$

$=\omega_n^p F[k+p] \sum_{i=0}^{n-1}\omega_n^i$

$=0$

所以 $nF[k]=\sum_{i=0}^{n-1}(\omega_n^{-k})^iG[i]$。

3.实现

void fft( Complex *F , int len , int op ) {
    if( len == 1 ) return;
    
    Complex *Fl = F , *Fr = F + len / 2;
    for( int i = 0 ; i < len ; i ++ ) tmp[ i ] = F[ i ];
    for( int i = 0 ; i < len / 2 ; i ++ )
        Fl[ i ] = tmp[ i << 1 ] , Fr[ i ] = tmp[ i << 1 | 1 ];
    
    fft( Fl , len / 2 , op );
    fft( Fr , len / 2 , op );

    Complex w = { cos( 2 * pi / len ) , op * sin( 2 * pi / len ) } , buf = { 1 , 0 };
    for( int i = 0 ; i < len / 2 ; i ++ ) {
        tmp[ i ] = Fl[ i ] + buf * Fr[ i ]; //*
        tmp[ i + len / 2 ] = Fl[ i ] - buf * Fr[ i ];
        buf = buf * w;
    }
    for( int i = 0 ; i < len ; i ++ ) F[ i ] = tmp[ i ];
}

很显然，递归版本因为常数太大 $T$ 掉了，考虑优化常数。

注意到 $(*)$ 处的 $buf*Fr[i]$ 计算了两次，但复数乘法极慢，可以用变量存下来。

但还是过不了，我们需要一种常数更小的写法。

五.快速傅里叶变换

看一下离散傅里叶变换中的系数变化：

$0,1,2,3,4,5,6,7$

$0,2,4,6,1,3,5,7$

$0,4,2,6,1,5,3,7$

把它们的二进制写出来，即

$\begin{cases} 0(000) \rightarrow 0(000) \\ 1(001) \rightarrow 4(100) \\ 2(010) \rightarrow 2(010) \\ 3(011) \rightarrow 6(110) \\ 4(100) \rightarrow 1(001) \\ 5(101) \rightarrow 5(101) \\ 6(110) \rightarrow 3(011) \\ 7(111) \rightarrow 7(111) \end{cases}$

经过观察发现，原数列位置的二进制反转后得到现在的位置。我不会证

递推公式也不很难想到：

$filp[i]=(filp[i>>1]>>1)|((i\&1)?limit>>1:0)$

然后就可以迭代实现了。

六.例题

1.P3803 【模板】多项式乘法

模板题，直接上代码（常数巨大）。

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;

const int MAXN = 2097152;
const double pi = acos( -1 );
struct Complex {
    double x , y;
	Complex(){}
	Complex( double X , double Y ) {
		x = X , y = Y;
	}

    Complex operator + ( const Complex &a ) const {
        return Complex( x + a.x , y + a.y );
    }
    Complex operator - ( const Complex &a ) const {
        return Complex( x - a.x , y - a.y );
    }
    Complex operator * ( const Complex &a ) const {
        return Complex( x * a.x - y * a.y , x * a.y + y * a.x );
    }
	Complex operator / ( const double &a ) const {
		return Complex( x / a , y / a );
	}
}f[ MAXN + 5 ] , g[ MAXN + 5 ] , h[ MAXN + 5 ] , tmp[ MAXN + 5 ];
int n , m , Filp[ MAXN + 5 ];

void fft( Complex *F , int op ) {
    for( int i = 0 ; i < n ; i ++ )
        if( i < Filp[ i ] ) swap( F[ i ] , F[ Filp[ i ] ] );

    for( int len = 2 ; len <= n ; len <<= 1 ) {
        Complex w( cos( 2 * pi / len ) , op * sin( 2 * pi / len ) );
        for( int l = 0 ; l < n ; l += len ) {
            Complex buf( 1 , 0 );
            for( int i = l ; i < l + len / 2 ; i ++ ) {
                Complex t = buf * F[ i + len / 2 ];
                F[ i + len / 2 ] = F[ i ] - t;
                F[ i ] = F[ i ] + t;
                buf = buf * w;
            }
        }
    }

	if( op == -1 )
		for( int i = 0 ; i <= n ; i ++ )
			F[ i ] = F[ i ] / n;
}

int main() {
    scanf("%d %d",&n,&m);
    for( int i = 0 ; i <= n ; i ++ ) scanf("%lf",&f[ i ].x);
    for( int i = 0 ; i <= m ; i ++ ) scanf("%lf",&g[ i ].x);
    for( m += n , n = 1 ; n <= m ; n <<= 1 );
    for( int i = 0 ; i < n ; i ++ )
        Filp[ i ] = ( Filp[ i >> 1 ] >> 1 ) | ( ( i & 1 ) ? n >> 1 : 0 );
    
    fft( f , 1 ) , fft( g , 1 );
    for( int i = 0 ; i < n ; i ++ ) h[ i ] = f[ i ] * g[ i ];
    fft( h , -1 );

    for( int i = 0 ; i <= m ; i ++ )
        printf("%d ", int( h[ i ].x + 0.49 ) );
    return 0;
}

2.P3338 [ZJOI2014]力

$E_j=\sum_{i=1}^{j-1} \frac{q_i}{(i-j)^2} - \sum_{i=j+1}^n \frac{q_i}{(i-j)^2}$

$E_j=\sum_{i=1}^{j} \frac{q_i}{(i-j)^2} - \sum_{i=j}^n \frac{q_i}{(i-j)^2}$

令 $f(i)=q_i , g(i)=(\frac{1}{i^2})$，特殊的，$f(0)=g(0)=0$

$E_j=\sum_{i=1}^j f(i) \times g(i-j) - \sum_{i=j}^n f(i) \times g(i-j)$

$E_j=\sum_{i=0}^j f(i) \times g(i-j) - \sum_{i=j}^n f(i) \times g(i-j)$

$E_j=\sum_{i=0}^j f(i) \times g(i-j) - \sum_{i=0}^{n-j} f(i+j) \times g(i)$

令 $f'[i]=f[n-i]$。

$E_j=\sum_{i=0}^j f(i) \times g(i-j) - \sum_{i=0}^{n-j} f(n-(i+j)) \times g(i)$

$E_j=\sum_{i=0}^j f(i) \times g(i-j) - \sum_{i=0}^{n-j} f((n-j)-i) \times g(i)$

直接计算卷积做差即可。

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int MAXN = 4e5; 
const double pi = acos( -1 );
struct Complex {
    double x , y;
	Complex(){}
	Complex( double X , double Y ) {
		x = X , y = Y;
	}

    Complex operator + ( const Complex &a ) const {
        return Complex( x + a.x , y + a.y );
    }
    Complex operator - ( const Complex &a ) const {
        return Complex( x - a.x , y - a.y );
    }
    Complex operator * ( const Complex &a ) const {
        return Complex( x * a.x - y * a.y , x * a.y + y * a.x );
    }
	Complex operator / ( const double &a ) const {
		return Complex( x / a , y / a );
	}
}f[ MAXN + 5 ] , rf[ MAXN + 5 ] , g[ MAXN + 5 ] , h[ MAXN + 5 ];
int n , l , Filp[ MAXN + 5 ];

void fft( Complex *F , int op , int n ) {
    for( int i = 0 ; i < n ; i ++ )
        if( i < Filp[ i ] ) swap( F[ i ] , F[ Filp[ i ] ] );
    for( int len = 2 ; len <= n ; len <<= 1 ) {
        Complex w( cos( 2 * pi / len ) , op * sin( 2 * pi / len ) );
        for( int l = 0 ; l < n ; l += len ) {
            Complex buf( 1 , 0 );
            for( int i = l ; i < l + len / 2 ; i ++ ) {
                Complex t = buf * F[ i + len / 2 ];
                F[ i + len / 2 ] = F[ i ] - t;
                F[ i ] = F[ i ] + t;
                buf = buf * w;
            }
        }
    }
	if( op == -1 )
		for( int i = 0 ; i < n ; i ++ )
			F[ i ] = F[ i ] / n;
}

int main() {
	scanf("%d",&n);
	for( int i = 1 ; i <= n ; i ++ ) {
		scanf("%lf", &f[ i ].x );
		rf[ n - i ] = f[ i ];
		g[ i ].x = 1.0 / i / i;
	}
	
	for( l = 1 ; l <= 2 * n ; l <<= 1 );
	for( int i = 0 ; i < l ; i ++ )
		Filp[ i ] = ( Filp[ i >> 1 ] >> 1 ) | ( ( i & 1 ) ? l >> 1 : 0 );
	
	fft( f , 1 , l ) , fft( g , 1 , l ) , fft( rf , 1 , l );
	for( int i = 0 ; i <= l ; i ++ ) f[ i ] = f[ i ] * g[ i ];
	for( int i = 0 ; i <= l ; i ++ ) rf[ i ] = rf[ i ] * g[ i ];
	fft( f , -1 , l ) , fft( rf , -1 , l );
	
	for( int i = 1 ; i <= n ; i ++ )
		printf("%.3f\n", f[ i ].x - rf[ n - i ].x );
    return 0;
}

3.P3763 [TJOI2017]DNA

$fft$ 求字符串匹配的套路，将 模式串 $T$ 反转，对 $A,T,C,G$ 分别求贡献。

要使相同位置字符不同的权值为一，将文本串中等于当前匹配字符的位置赋为 $1$ ，模式串中不等于当前字符的赋为 $1$。

最后卷积和权值小于三的便是匹配成功的位置。

这道题 $fft$ 常数过大需要吸氧。

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int MAXN = 4e5; 
const double pi = acos( -1 );
struct Complex {
    double x , y;
	Complex(){}
	Complex( double X , double Y ) {
		x = X , y = Y;
	}

    Complex operator + ( const Complex &a ) const {
        return Complex( x + a.x , y + a.y );
    }
    Complex operator - ( const Complex &a ) const {
        return Complex( x - a.x , y - a.y );
    }
    Complex operator * ( const Complex &a ) const {
        return Complex( x * a.x - y * a.y , x * a.y + y * a.x );
    }
	Complex operator / ( const double &a ) const {
		return Complex( x / a , y / a );
	}
}f[ MAXN + 5 ] , g[ MAXN + 5 ] , h[ MAXN + 5 ];
char Hash[ 4 ] = { 'A' , 'T' , 'C' , 'G' };
int t , n , m , l , p[ MAXN + 5 ] , Filp[ MAXN + 5 ];
string S , T;

void print( Complex *F ) {
	for( int i = 0 ; i < l ; i ++ )
		printf("%lf %lf\n",F[i].x,F[i].y);
	printf("\n");
}
void fft( Complex *F , int op , int n ) {
	//print( F );
    for( int i = 0 ; i < n ; i ++ )
        if( i < Filp[ i ] ) swap( F[ i ] , F[ Filp[ i ] ] );

	//print( F );

    for( int len = 2 ; len <= n ; len <<= 1 ) {
        Complex w( cos( 2 * pi / len ) , op * sin( 2 * pi / len ) );
        for( int l = 0 ; l < n ; l += len ) {
            Complex buf( 1 , 0 );
            for( int i = l ; i < l + len / 2 ; i ++ ) {
                Complex t = buf * F[ i + len / 2 ];
                F[ i + len / 2 ] = F[ i ] - t;
                F[ i ] = F[ i ] + t;
                buf = buf * w;
            }
        }
    }

	//print( F );

	if( op == -1 )
		for( int i = 0 ; i < n ; i ++ )
			F[ i ] = F[ i ] / n;
}

int main() {
	scanf("%d",&t);
	while( t -- ) {
		memset( p , 0 , sizeof( p ) );

		cin >> S >> T;
		reverse( T.begin( ) , T.end( ) );
		n = S.length( ) , m = T.length( ) , l = 1;
		
		for( ; l < n + m ; l <<= 1 );
		for( int i = 0 ; i < l ; i ++ )
			Filp[ i ] = ( Filp[ i >> 1 ] >> 1 ) | ( ( i & 1 ) ? l >> 1 : 0 );
		
		for( int c = 0 ; c < 4 ; c ++ ) {
			for( int i = 0 ; i < l ; i ++ ) f[ i ] = Complex( i < n ? ( S[ i ] != Hash[ c ] ) : 0 , 0 );
			for( int i = 0 ; i < l ; i ++ ) g[ i ] = Complex( i < m ? ( T[ i ] == Hash[ c ] ) : 0 , 0 );
			//print( f );
			fft( f , 1 , l ) , fft( g , 1 , l );
			for( int i = 0 ; i < l ; i ++ ) h[ i ] = f[ i ] * g[ i ];
			fft( h , -1 , l );
			
			for( int i = m - 1 ; i < n ; i ++ )
				p[ i ] += int( h[ i ].x + 0.5 );
		}

		int Ans = 0;
		for( int i = m - 1 ; i < n ; i ++ )
			Ans += ( p[ i ] <= 3 );
		printf("%d\n", Ans );
	}
    return 0;
}

参考资料

洛谷日报